unit tile_mc;

interface
uses common;
Var
   nodum : boolean;
procedure init_th(Var epot :myreal);
procedure mc_tiling(swo,cpl,vpin,lvl0,tlvl,kk,mcs,lab: integer; Var tpot,rate:myreal);

implementation
uses d04_header,ppot,tools,etools,outp,decor;
Type
   b50 = array[1..50] of boolean;
   i02 = array[0..2] of integer;
Var
   invers,mirror,identic : array[1..4] of boolean;
   txerr		 : text;
Const
   atlist = 1000;
   {for local swap-pair}
   iucp	 = 4;
   iucpp = (iucp*2+1)*(iucp*2+1)*(iucp*2+1); {for flat lists gen. by 3 iext loops}
Var
   {these arrays are buffers for efficient energy update.}
   r_0,r_1,r_2 : array[1..atlist] of v3;
   t_0,t_1,t_2 : array[1..atlist] of integer;
   po_1        : array[1..atlist] of i2;
   n_0,n_1,n_2 : integer;
   {for local swap-pair}
   {tt1,tt2     : array[1..atlist,1..atlist,1..t_ele,1..t_ele] of myreal;}
   npx         : integer;
   px          : array[1..iucpp] of v3;

procedure table_update(ii,nnl:integer; fln:flist);
Var
   i,j,k,kk,o,io  : integer;
   j1,j2,ds,i1,i2 : integer;
   d2,eps,xx	  : myreal;
   u,v,w	  : v2;
begin
   d2:=1/tau/tau; {node-node distance for sigmaflip}
   eps:=0.001; {accuracy to match bonds}
   i1:=0; i2:=0;
   {writeln('tab-upd ',ii:1,' ',nnl:1);}
   if (tnodes0<50) then begin i1:=1; i2:=1; end; {for small unit cells}
   if((io3d>0) and (multi>1)) then begin
      for io:=-1 to 1 do begin
	 for o:=0 to 10 do ltm[ii,io,o]:=0;
	 nvm[ii,io]:=0;
	 ltc[ii,io]:=0;
	 for k:=1 to nnl do begin
	    j:=fln[k];
	    for o:=0 to 10 do if (ltm[j,io,o]=ii) then begin
	       ltm[j,io,o]:=0;
	       if(io>0) then nvm[j,io]:=nvm[j,io]-1;
	    end;
	    if(ltc[j,io]=ii) then ltc[j,io]:=0;
	 end;
      end;
   end;
   nv[ii]:=0;
   for j:=1 to 10 do lt[ii,j]:=0;
   for k:=1 to nnl do begin
      j:=fln[k];
      for o:=1 to 10 do if (lt[j,o]=ii) then begin
	 nv[j]:=nv[j]-1;
	 lt[j,o]:=0;
      end;
   end;
   for kk:=1 to nnl do begin
      j:=fln[kk];
      ds:=t2t(ii,j);
      if abs(ds)<2 then begin
	 for j1:=-i1 to i1 do
	    for j2:=-i2 to i2 do begin
	       for o:=1 to 2 do v[o]:=bound(rt[j,o]-rt[ii,o]);
	       v[1]:=v[1]+j1; v[2]:=v[2]+j2;
	       ortho(v,u);
	       if (absv(u)<eps) then begin
		  ltm[ii,ds,0]:=j;
		  ltm[j,-ds,0]:=ii;
		  if (multi=2) then begin {periodic boundary for multi=2}
		     ltm[ii,-ds,0]:=j;
		     ltm[j,ds,0]:=ii;
		  end;
	       end;
	       {this is for sigma-flip in two-level tiling}
	       if sigmaflip then if (abs(absv(u)-d2)<eps) then begin
		  ltc[ii,ds]:=j;
		  ltc[j,-ds]:=ii;
		  if (multi=2) then begin
		     ltc[ii,-ds]:=j;
		     ltc[j,ds]:=ii;
		  end;
	       end;
	       if (abs(absv(u)-1.0)<eps) then begin
		  for k:=1 to 10 do begin
		     xx:=sp(u,re[k]);
		     if (abs(xx-1.0)<eps) then begin
			nvm[ii,ds]:=nvm[ii,ds]+1;
			ltm[ii,ds,k]:=j;
			nvm[j,-ds]:=nvm[j,-ds]+1;
			ltm[j,-ds,po(k+5)]:=ii;
			if (multi=2) then begin {periodic boundary for multi=2}
			   nvm[ii,-ds]:=nvm[ii,-ds]+1;
			   ltm[ii,-ds,k]:=j;
			   nvm[j,ds]:=nvm[j,ds]+1;
			   ltm[j,ds,po(k+5)]:=ii;
			end;
			if(ds=0) then begin
			   lt[ii,k]:=j;
			   nv[ii]:=nv[ii]+1;
			   lt[j,po(k+5)]:=ii;
			   nv[j]:=nv[j]+1;
			end;
		     end;
		  end;
	       end;
	    end;
      end;
   end;
   if((io3d>0) and (multi>1)) then begin
      for o:=-1 to 1 do if (o<>0) then begin
	 if((ltc[ii,o]=0) and (ltm[ii,o,0]=0)) then begin
	    writeln('table_update ',ii:1,' : LTC error');
	    halt
	 end;
      end;
   end;
end; { table_update }

procedure init_th(Var epot : myreal);
var j,k,l,o : integer;
begin
   epot:=0.0;
   for o:=1 to n_obj do begin
      obj_st[o]:=0;
      for j:=1 to t_sym do obj_st_list[o,j]:=0;
   end;
   for o:=1 to n_objch do obj_sta[o]:=0;
   for o:=1 to n_objch do objcs[o]:='';
   for o:=1 to n_obj do if obj_code[o]<99 then
      obj_stn[obj_code[o]]:=round(obj_stx[obj_code[o]]*tnodes); {last obj_stx of given obj_code counts}
   if(instances) then for l:=1 to tnodes do new_pat_instance(l,0);
   for l:=1 to tnodes do new_node_objects(l,0);
   for l:=1 to tnodes do begin
      for j:=1 to n_objmu[l] do begin
	 o:=n_objpo[l,j,1]; {object}
	 k:=n_objpo[l,j,2]; {orientation}
	 obj_st[o]:=obj_st[o]+1;
	 obj_st_list[o,k]:=obj_st_list[o,k]+1;
	 if (obj_code[o]<99) then obj_sta[obj_code[o]]:=obj_sta[obj_code[o]]+1;
	 epot:=epot+obj_ham[o];
      end;
   end;
   for o:=1 to n_obj do if (obj_code[o]<99) then begin
      objcs[obj_code[o]]:=objcs[obj_code[o]]+obj_name[o];
   end;
   for o:=1 to n_objch do writeln('Constrained obj-freq mode(curr,req,chp): ',obj_sta[o]:1,' ',obj_stn[o]:1,' ',obj_ch[o]:1:3,' ',objcs[o]);
end; { init_th }
procedure E_update_thru_nodes(mu1,mu : integer;x:flist; Var u0,u1:myreal);
var i,j,k : integer;
   u	  : myreal;
begin
   {This is standard and slow algorithm of calculating dE of the
    configurations before and after flip.
    Atomic flips (reflected in node.atom constructs) occure for
    7 nodes (flippable + 6 neighbors of hexagon).
    We need to consider more nodes within Rcut. All _mu_ nodes
    needing update are in the list _x_.
    The old record of the 7 nodes is stored using pointer TNODES,
    behind the last actual record. We need to compute:
    1. "self" energies of the first 7 nodes, for old/new configs.
    2. node-node energies for the first 7 nodes
    3. node-node energies between the first 7 and other nodes in
    the list _x_. The following ocmputes all that in a compact form.}
   u1:=0.0; u0:=0.0;
   for i:=1 to mu1{7} do
      for j:=i to mu do begin
	 node_pair_energy(x[i],x[j],u); u1:=u1+u;
	 {there are 7 old nodes stored on the stack "tnodes+i".}
	 if (j<mu1+1) then k:=tnodes+j else k:=x[j];
	 node_pair_energy(tnodes+i,k,u); u0:=u0+u;
      end;
end; { update_through_nodes }
procedure flip_union(eps : myreal);
var
   i,j,o       : integer;
   x	       : v3;
   por	       : myreal;
   b_0,b_1,b_2 : array[1..atlist] of boolean;
begin
   if (debug[3]) then write('out of ',n_1:1,'/',n_0:1,' atoms, ');
   {for i:=1 to n_0 do begin
   if (t_0[i]<>t_1[i]) then writeln('FU: chem mismatch');
   write(i:1,' ',t_0[i]:1,'/',t_1[i]:1,'  ');
   for o:=1 to 3 do write(r_0[i,o]:1:3,' '); write('  ');
   for o:=1 to 3 do write(r_1[i,o]:1:3,' '); writeln;
   end;}
   for i:=1 to n_0 do b_0[i]:=true;
   for i:=1 to n_1 do b_1[i]:=true;
   {identify atoms that actually stay fixed during the flip}
   for i:=1 to n_1 do if b_1[i] then
      for j:=1 to n_0 do if b_0[j] then if (t_1[i]=t_0[j]) then begin
	 por:=dist3(r_1[i],r_0[j]);
	 if (por<eps) then begin
	    b_1[i]:=false; b_0[j]:=false;
	    {calculate mean of these two}
	    {for o:=1 to 3 do x[o]:=(bound(r_0[j,o])+bound(r_1[i,o]))/2.0;
	    r_1[i]:=x; r_0[j]:=x;}
	 end;
      end;
   {move these atoms to the end of the r_2 list, we do this for
    BOTH r_0 and r_1 to preserve symmetry, when there are very
    small differences in positions}
   j:=0; for i:=1 to n_1 do if b_1[i] then begin
      j:=j+1; r_1[j]:=r_1[i]; t_1[j]:=t_1[i];
   end else begin n_2:=n_2+1;
      r_2[n_2]:=r_1[i]; t_2[n_2]:=t_1[i];
   end;
   n_1:=j;
   j:=0; for i:=1 to n_0 do if b_0[i] then begin
      j:=j+1; r_0[j]:=r_0[i]; t_0[j]:=t_0[i];
   end;
   { else begin n_2:=n_2+1;
	r_2[n_2]:=r_0[i]; t_2[n_2]:=t_0[i];
   end;}
   n_0:=j;
   if debug[3] then begin
      writeln(n_1:1,'/',n_0:1,' actually flipped; ',n_2:1,' fixed atoms.');
   end;
end;
procedure add2list(po:integer; Var mu:integer; Var x:flist);
var i,j,k,k0,o,n2 : integer;
   ar,ar0	  : myreal;
   found	  : boolean;
   rc		  : myreal;
begin
   rc:=rcutn; {rcutn is common from d04_header}
   for i:=1 to tnodes do if (i<>po) then begin
      ar:=dist3(rt[po],rt[i]);
      if (ar<rc) then begin
	 found:=false; for j:=1 to mu do if (x[j]=i) then found:=true;
	 if not(found) then begin
	    mu:=mu+1;
	    if (mu>t_flip)then begin
	       writeln('add2list: Insufficient array dimension t_flip ',t_flip:1,' in d04_header.p.');
	       halt;
	    end;
	    x[mu]:=i;
	 end;
      end;
   end;
end; { add2list }
procedure add_boundary(dep : integer; Var mu:integer; Var x:flist);
{Note: this only adds 2D boundary vertices!!}
Var
   i,ll,j,o,n0,n1 : integer;
   b		  : array[1..tdim] of boolean;
begin
   for i:=1 to tnodes do b[i]:=false;
   for i:=1 to mu do b[x[i]]:=true;
   n0:=1; n1:=mu;
   for i:=1 to dep do begin
      n1:=mu;
      for j:=n0 to n1 do begin
	 for o:=1 to 10 do begin
	    ll:=lt[x[j],o];
	    if (ll>0) and not(b[ll]) then begin
	       mu:=mu+1; x[mu]:=ll; b[ll]:=true;
	    end;
	 end;
      end;
      n0:=n1+1;
   end;
end; { add_boundary }
procedure E_update_buf(mu0,mu: integer;x:flist; Var u0,u1:myreal);
var i,j	: integer;
   u	: myreal;
   ide	: boolean;
begin
   {This is fast algorithm of calculating dE before and after flip.
    r_0 is list of deleted atoms.
    r_1 is list of new atoms.
    r_2 is list of fixed atoms.
    [these are stored inside _flip()_ when buf=true], and are
    declared in flip.p .
    We need to update/compute:
    1. E of r_0 and r_1 interacting with r_2
    2. E of r_0 and r_1 interacting with atoms bound to Flist of nodes,
    starting from position 8 (flippable vertex=1, hexagon 2-7)
    3. "self" Energies of lists r_0 and r_1.
    Another speedup: some new/old atoms are actually identical.
    These may be moved from r_0/r_1 to r_2.
    NOTE: unlike node.atoms, atom lists r_x already HAVE the
    node-coordinate in them.}
   {compute the energies u0 and u1}
   u0:=0.0; u1:=0.0;
   {"self" energies}
   for i:=1 to n_0 do
      for j:=i to n_0 do begin
	 ide:=(i=j);
	 atom_pair_energy(ide,r_0[i],r_0[j],t_0[i],t_0[j],u);
	 {writeln('self ',i:1,'-',j:1,' ',t_0[i],' ',t_0[j],' ',u:1:4);}
	 u0:=u0+u;
	 if (i=j) then u0:=u0-u/2.0;
      end;
   for i:=1 to n_1 do
      for j:=i to n_1 do begin
	 ide:=(i=j);
	 atom_pair_energy(ide,r_1[i],r_1[j],t_1[i],t_1[j],u);
	 u1:=u1+u;
	 if (i=j) then u1:=u1-u/2.0;
      end;
   {writeln('self u1 u0 ',u1:1:4,' ',u0:1:4);}
   {pair energy with the list r_2}
   ide:=false;
   for i:=1 to n_2 do begin
      for j:=1 to n_0 do begin
	 atom_pair_energy(ide,r_2[i],r_0[j],t_2[i],t_0[j],u);
	 u0:=u0+u;
      end;
      for j:=1 to n_1 do begin
	 atom_pair_energy(ide,r_2[i],r_1[j],t_2[i],t_1[j],u);
	 u1:=u1+u;
      end;
   end;
   {writeln('self+pair u1 u0 ',u1:1:4,' ',u0:1:4);}
   {finally, we need the atom-node energies}
   for i:=mu0+1 to mu do begin
      for j:=1 to n_0 do begin
	 atom_node_energy(x[i],t_0[j],r_0[j],u);
	 u0:=u0+u;
      end;
      for j:=1 to n_1 do begin
	 atom_node_energy(x[i],t_1[j],r_1[j],u);
	 u1:=u1+u;
      end;
   end;
   {writeln('final u1 u0 ',u1:1:4,' ',u0:1:4);}
end; { E_update_buf }
procedure export_replay(jj,ii:integer);
Var k,o,l,m : integer;
   noexport   : boolean;
begin
   m:=0; for o:=1 to n_obj do
      m:=m+obj_st[o]+obj_st[o+k]-obj_st_u[o]-obj_st_u[o+k];
   noexport:=(nodum and (abs(dfene[jj])<0.0001) and (m=0));
   {writeln(jj:1,' ',m:1,' ',abs(fene[jj]):1:5,' ',noexport);}
   noexport:=(noexport or ((fene[jj]<erec1) and (fene[jj]>erec2)));
   {noexport:=(doflip[jj]=0);} {quick-fix!}
   if not(noexport) then begin
      write(txfout,jj:5,' ',ii:4,' ');
      for o:=1 to n_obj do begin
	 if dfe then begin
	    l:=obj_st[o]-obj_st_u[o]
	 end else begin
	    l:=obj_st[o];
	 end;
	 write(txfout,' ',l:2);
      end;
      if dfe then writeln(txfout,'   ',dfene[jj]:1:6) else
	 writeln(txfout,'   ',fene[jj]:1:6);
   end;
end; { export_replay }
procedure copy_nrec(ifrom,ito :integer);
begin
   rt[ito]:=rt[ifrom];
   bt[ito]:=bt[ifrom];
   nu[ito]:=nu[ifrom];
   {tvl[ito]:=tvl[ifrom];
   tvlpo[ito]:=tvlpo[ifrom];}
end;
procedure copy_trec(ifrom,ito:integer);
var i,o,io : integer;
begin
   nv[ito]:=nv[ifrom];
   lt[ito]:=lt[ifrom];
   {if (multi>1) then io:=1 else io:=0;}
   for o:=-io3d to io3d do begin
      nvm[ito,o]:=nvm[ifrom,o];
      ltm[ito,o]:=ltm[ifrom,o];
      ltc[ito,o]:=ltc[ifrom,o];
   end;
   for i:=1 to n_objmu[ifrom] do begin
      n_objpo[ito,i]:=n_objpo[ifrom,i];
      n_obj_fl[ito,i]:=n_obj_fl[ifrom,i];
      n_obj_cl[ito,i]:=n_obj_cl[ifrom,i];
   end;
   n_objmu[ito]:=n_objmu[ifrom];
   for o:=1 to tdim do n2n[ifrom,o]:=n2n[ito,o];
end; { copy_trec }
procedure copy_arec(ifrom,ito : integer);
var i,o : integer;
begin
   for i:=1 to nan[ifrom] do begin
      n_r[ito,i]:=n_r[ifrom,i];
      if repre5 then n_r5[ito,i]:=n_r5[ifrom,i];
      n_para[ito,i]:=n_para[ifrom,i];
      n_a[ito,i]:=n_a[ifrom,i];
   end;
   nan[ito]:=nan[ifrom];
end; { copy_arec }
procedure restore_rec(nfl,nft: integer; fln,flt:flist);
Var
   i,j,k : integer;
begin
   if debug[2] then write('restore_rec (',nfl:1,')  ');
   for i:=1 to nfl do begin
      if (debug[2]) then begin
	 k:=tvl[fln[i]];
	 write(' R',fln[i]:1,' ');
      end;
      copy_nrec(tnodes+i,fln[i]);
      copy_trec(tnodes+i,fln[i]);
      copy_arec(tnodes+i,fln[i]);
      if (debug[2]) then begin
      end;
   end;
   for i:=1 to nft do begin
      j:=i+nfl;
      if (debug[2]) then begin
	 write(' A',flt[i]:1,':',tvl[flt[i]]:1,'>');
      end;
      copy_nrec(tnodes+j,flt[i]);
      copy_trec(tnodes+j,flt[i]);
   end;
   if debug[2] then writeln;
end; { restore_rec }
procedure fnodelist(Var mu : integer; Var x:flist; po:integer);
var o,i	 : integer;
   found : boolean;
begin
   {writeln(dp:1,' ',po:1,' ',mu:1,' ');}
   for o:=1 to 10 do if (lt[po,o]>0) then begin
      found:=false; for i:=1 to mu do
	 if not(found) then if (x[i]=lt[po,o]) then found:=true;
      if not(found) then begin
	 mu:=mu+1; x[mu]:=lt[po,o];
	 if (mu=t_flip) then begin
	    writeln('Insufficient array dim for t_flip=',t_flip:1);
	    halt;
	 end;
      end;
   end;
end; { fnodelist }
procedure nodes_union(Var mu: integer; Var x:flist);
var i,j,mu0 : integer;
   x0	    : flist;
   add	    : boolean;
begin
   mu0:=1; x0[1]:=x[1];
   for i:=2 to mu do begin
      add:=true;
      for j:=1 to mu0 do if (add and (x[i]=x0[j])) then add:=false;
      if add then begin mu0:=mu0+1; x0[mu0]:=x[i]; end;
   end;
   mu:=mu0; x:=x0;
end;
procedure flipit(jfl:integer;flst:flist ;Var nnmu:integer;Var fln:flist;Var done:boolean);
Var
   i,j,k,l,o,ii,nux,dep	: integer;
   p1,p2,nu0,nu1,kl	: integer;
   nsurf,nnin,nini	: integer;
   nstart,nend,nsw,nuii	: integer;
   refdir,refpo,dir	: integer;
   sg,sgc,fc,nuc	: integer;
   refdirx		: integer;
   level,depth		: integer;
   dn			: array[1..6] of i2;
   dnm			: array[1..6,-1..1] of i2;
   ijk			: array[1..2,2..3] of integer;
   fv			: i5;
   path			: vpath;
   ix,lz		: array[-1..1] of integer;
   zflippable		: boolean;
   ident,invt,dbx	: boolean;
   w			: v3;
   fff			: nodelist;
procedure fnode(var dp	: integer; p:vpath; var l:integer);
var i: integer;
begin
   i:=l;
   l:=lt[i,p[dp]];
   dp:=dp-1;
   if (dp>0) then fnode(dp,p,l);
end; { fnode }
begin
   ii:=flst[jfl];
   nini:=nnmu;
   nnmu:=nnmu+1; fln[nnmu]:=ii; nnin:=nnmu;
   nsurf:=6;
   done:=false;
   ijk[1,2]:=3;
   ijk[1,3]:=2;
   ijk[2,3]:=1;
   {3 old neighbors and 3 old directions}
   k:=0; for j:=1 to 10 do if (lt[ii,j]>0) then begin
      k:=k+1;
      dn[k,1]:=lt[ii,j]; {neighbor}
      dn[k,2]:=j; {direction}
   end;
   if (k<>3) then begin
      write('FLIPIT (',ii:1,'): NV<>3 (stated ',nv[ii]:1,'): ');
      for i:=1 to 10 do write(lt[ii,i]:1,' '); writeln;
      if (flst[0]>1) then writeln('flst: ',flst[1]:1,' ',flst[2]:1);
      writeln('Exporting config into copa.err.');
      export_config_fl('copa.err');
      halt;
   end;
   {The following procedure implicitly finds inv-symmetry related
    nodes, such that 1-4 2-5 3-6 are related. This property will
    be used to swap the site list}
   refdir:=0;
   for i:=1 to 2 do {3 new neighbors and directions}
      for j:=i+1 to 3 do begin
	 path[1]:=dn[i,2]; path[2]:=dn[j,2];
	 l:=ii; {we start from THIS node.}
	 dep:=2; {depth of the path}
	 fnode(dep,path,l); {finds node label through vector path}
	 k:=ijk[i,j]; {third node/direction out of the 3}
	 dn[k+3,1]:=l;
	 dn[k+3,2]:=dn[k,2];
	 p1:=po(dn[i,2]-dn[k,2]); if (p1>5) then p1:=p1-10;
	 p2:=po(dn[j,2]-dn[k,2]); if (p2>5) then p2:=p2-10;
	 if (abs(p1)=abs(p2)) then begin
	    refdir:=dn[k,2];
	    refpo:=dn[k,1];
	 end;
      end; { fnode }
   for o:=1 to nsurf do fln[nnin+o]:=dn[o,1];
   nnmu:=nnin+nsurf;
   zflippable:=true;
   if (io3d>0) then begin
      l:=nnmu; {7; this tiling nodelist already set}
      for o:=-1 to 1 do lz[o]:=0; 
      for o:=-1 to 1 do if (o<>0) then begin
	 if (node_z_alligned(ii,o)<>ltm[ii,o,0]) then begin
	    writeln('flipit LTM problem: ',ltm[ii,o,0]:1,' ',node_z_alligned(ii,o):1);
	    halt;
	 end;
	 lz[o]:=ltm[ii,o,0]; {now implemented directly into the tables}
	 if (lz[o]=0) then begin {the top/bottom non-z-alligned vertex}
	    if (ltc[ii,o]=0) then begin
	       writeln('Z-connectivity error ',ii:1,': ltm=0 AND ltc=0.');
	       halt;
	    end;
	    l:=l+1; fln[l]:=ltc[ii,o];
	 end else begin
	    l:=l+1; fln[l]:=lz[o];
	 end;
	 for k:=1 to nsurf do begin
	    kl:=k;
	    if (l=0) then begin
	       if (k<4) then kl:=k+3;
	       if (k>3) then kl:=k-3;
	    end;
	    dnm[kl,o,1]:=ltm[dn[k,1],o,0]; {pointer to node}
	    dnm[kl,o,2]:=dn[k,2]; {direction}
	    l:=l+1; fln[l]:=dnm[kl,o,1];
	    if (dnm[kl,o,1]=0) then zflippable:=false;
	 end;
      end;
      if (flip_man and not(zflippable)) then begin
	 writeln('Vertex ',ii:1,' is not z-flippable.');
      end;
      {if ((multi>1) and (l<>21)) then begin
      writeln('Z-stuff: wrong number of neighbors (',l:1,' for update (<>21)');
      halt;
      end;}
      nnmu:=l;
   end; { (multi>1) stuff }
   {5D flip vector}
   fv:=bt[ii];
   for o:=1 to 5 do for j:=1 to 3 do fv[o]:=fv[o]+bs[dn[j,2],o];
   bound5(fv,w);
   nux:=dperp_offset(-nu_offs,nui(fv));
   nuii:=dperp_offset(-nu_offs,nu[ii]);
   if (debug[2]) then begin
      writeln('-------------- flipit ---------------------------------');
      write('FLIPPING: ',ii:1,',',jfl:1,'/',flst[0]:1,' dperp[',nuii:1,'->',nux:1,']; flip_perp_allowed : ',flip_b[nuii,nux]);
      if (multi>1) then write(' z-flippable: ',zflippable);
      writeln;
      {write(flip_b[nu[ii],nux],' ',zflippable,' + ');
      writeln(nu[ii]:1,' ',nuii:1,' ',nux:1,' ',flip_b[nuii,nux],' ',(multi>1));}
   end;
   if (flip_b[nuii,nux] and zflippable) then begin {FLIP}
      if(instances) then begin
	 for o:=1 to tnodes do if n2n[o,flst[jfl]] then begin
	    if debug[2] then writeln('flipping vertex ',flst[jfl]:1,' forces update for vertex ',o:1);
	    {write('yyYY ');
	    for k:=1 to tdim do if n2n[o,k] then write(' ',k:1); writeln;}
	    new_pat_instance(o,1);
	    if(svn>0) then begin
	       if debug[2] then writeln('  ...flipit() ipat-dep ',o:1,' ',svn:1);
	       for j:=1 to svn do begin
		  for k:=1 to 10 do if (lt[svlist[j],k]>0) then begin
		     {writeln('XX ',svlist[j]:1,' ',lt[svlist[j],k]);}
		     nnmu:=nnmu+1; fln[nnmu]:=lt[svlist[j],k];
		  end;
	       end;
	    end;
	 end;
	 {fff[0]:=0;
	 for i:=nnin to nnmu do pat_instance_dep(fln[i],fff);
	 for i:=1 to fff[0] do begin nnmu:=nnmu+1; fln[nnmu]:=fff[i]; end;}
      end;
      nodes_union(nnmu,fln); {needed for combined flip, but also multi=2}
      if(debug[2]) then begin
	 write('  ...flipit() nnmu fln : ',nnmu:1,' ');
	 for j:=1 to nnmu do write(' ',fln[j]:1);
	 writeln;
      end;
      if deco then for i:=nnin to nnmu do COPY_AREC(fln[i],tnodes+i);
      for i:=nnin to nnmu do copy_nrec(fln[i],tnodes+i);
      for i:=nnin to nnmu do copy_trec(fln[i],tnodes+i);
      {optionally clear old pattern instances}
      if(instances) then begin {we do this BEFORE flip, for more vertices}
	 {for i:=nnin to nnmu do new_pat_instance(fln[i],1);}
      end;
      {write(flip_b[nuii,nux],' ',zflippable,'  -> ');}
      {nuc:=dperp_offset(-nu_offs,nu[ii]);}
      nuc:=nuii;
      if (debug[7]) then begin
	 write('LT(bef) ',ii:1,' -> ');
	 for o:=1 to 10 do if (lt[ii,o]>0) then write(lt[ii,o]:1,' '); writeln;
      end;
      {UPDATE ALGORITHM. We're updating 3-fold PT vertex.
       The flip vector is Sum_3_edges. Each of the 3 new bonds
       will be emanating from vertex departed by a sum of a pair of
       edges - from the old position of the 3-fold vertex -
       such that the bond itself is translated 3rd "old"
       bond}
      bt[ii]:=fv; {new vertex position}
      nu[ii]:=nui(fv);
      nuii:=dperp_offset(-nu_offs,nu[ii]);
      rt[ii]:=w;
      rt[ii,3]:=((ii-1) div tnodes0)/multi;
      for j:=1 to 3 do begin {erase old bonds}
	 p1:=dn[j,1]; p2:=dn[j,2];
	 if (lt[p1,po(p2+5)]<>ii) then {check} begin
	    write(ii:1,' : '); for o:=1 to 10 do write(lt[ii,o]:1,' '); writeln;
	    write(p1:1,' ',po(p2+5):1,' : ');
	    for o:=1 to 10 do write(lt[p1,o]:1,' '); writeln;
	    writeln('Erase-bonds: inconsistent update ',ii:1,'<>',lt[p1,po(p2+5)]:1); halt;
	 end;
	 lt[ii,p2]:=0;
	 ltm[ii,0,p2]:=0;
	 nv[ii]:=nv[ii]-1;
	 nvm[ii,0]:=nvm[ii,0]-1;
	 lt[p1,po(p2+5)]:=0;
	 ltm[p1,0,po(p2+5)]:=0;
	 nv[p1]:=nv[p1]-1;
	 nvm[p1,0]:=nvm[p1,0]-1;
      end;
      for k:=1 to 3 do begin {create new bonds}
	 l:=dn[k+3,1];
	 p2:=dn[k+3,2];
	 lt[l,p2]:=ii;
	 ltm[l,0,p2]:=ii;
	 nv[l]:=nv[l]+1;
	 nvm[l,0]:=nvm[l,0]+1;
	 lt[ii,po(p2+5)]:=l;
	 ltm[ii,0,po(p2+5)]:=l;
	 nv[ii]:=nv[ii]+1;
	 nvm[ii,0]:=nvm[ii,0]+1;
      end;
      if (nv[ii]<>3) then writeln('Update error: nv<>3.');
      if (debug[7]) then begin
	 write('LT (aft) ',ii:1,' -> ');
	 for o:=1 to 10 do if (lt[ii,o]>0) then write(lt[ii,o]:1,' '); writeln;
      end;
      if(instances) then begin
	 {for i:=nnin to nnmu do new_pat_instance(fln[i],0);}
      end;
      if (io3d>0) then begin 
	 {write('BEF ',lz[-1]:1,' ',lz[1]:1);}
	 for o:=-1 to 1 do if (o<>0) then begin
	    for j:=1 to 3 do begin {erase old bonds}
	       p1:=dnm[j,o,1]; p2:=dnm[j,o,2];
	       if (ltm[p1,-o,po(p2+5)]<>ii) then
		  writeln('Flipit: vertical connectivity problem: ',ii:1,' (',dperp_offset(-nu_offs,nu[p1]):1,') <> ',ltm[p1,-o,po(p2+5)]:1);
	       ltm[ii,o,p2]:=0;
	       nvm[ii,o]:=nvm[ii,o]-1;
	       ltm[p1,-o,po(p2+5)]:=0;
	       nvm[p1,-o]:=nvm[p1,-o]-1;
	    end;
	    for k:=1 to 3 do begin {create new bonds}
	       l:=dnm[k+3,o,1];
	       p2:=dnm[k+3,o,2];
	       ltm[l,-o,p2]:=ii;
	       nvm[l,-o]:=nvm[l,-o]+1;
	       ltm[ii,o,po(p2+5)]:=l;
	       nvm[ii,o]:=nvm[ii,o]+1;
	    end;
	    p1:=ltm[ii,o,0];
	    if (p1>0) then begin {ii and p1 were alligned: erase}
	       ltm[ii,o,0]:=0;
	       ltm[p1,-o,0]:=0;
	       ltc[ii,o]:=p1;
	       ltc[p1,-o]:=ii;
	    end else begin {ii and p1 NOT alligned: create z-allignement entry}
	       {p1:=node_z_alligned(ii,o);}
	       p1:=ltc[ii,o];
	       ltm[ii,o,0]:=p1;
	       ltm[p1,-o,0]:=ii;
	       ltc[ii,o]:=0;
	       ltc[p1,-o]:=0;
	    end;
	 end;
	 for o:=-1 to 1 do if (o<>0) then begin
	    lz[o]:=ltm[ii,o,0]; {update node z-allignments}
	 end;
	 {writeln('   AFT ',lz[-1]:1,' ',lz[1]:1);}
      end;
      done:=true;
   end; {flippable}
   if debug[2] then writeln('-------------- end_of_flipit -------------------');
end; {flipit}
procedure CMPCHL(swo :integer; flst :flist; nfl : integer;fln:flist );
{
"compare chemical lists":
find atoms that didnt change position after the flip (r_2) and need not
energy update); find pairs of atoms that do need E-update (r_0, r_1)
}
Var
   i0,i,j,k,o : integer;
   nm,nre,icp : integer;
   b,bx	      : array[1..at_of_node] of boolean;
   bfl	      : array[1..npth] of boolean;
begin
   n_0:=0; n_1:=0; n_2:=0;
   {set bfl true when node flipped}
   for i:=1 to nfl do begin
      {writeln('cmpchl-gagaga ',i:2,' ',nfl:1);}
      {if(swo=1) then bfl[i]:=false else begin}
	 bfl[i]:=false;
	 for j:=1 to flst[0] do if fln[i]=flst[j] then bfl[i]:=true;
	 {bfl[i]:=((fln[i]=flst[1]) or (fln[i]=flst[2]));}
      {end;}
   end;
   nre:=0;
   for i0:=1 to nfl do if not(bfl[i0]) then begin
      i:=fln[i0];
      icp:=tnodes+i0;
      nm:=0;
      for j:=1 to nan[i] do b[j]:=false;
      for j:=1 to nan[icp] do bx[j]:=false;
      {shouldnt do "flat" comparison first?}
      for j:=1 to nan[i] do begin
	 for k:=1 to nan[icp] do if not(bx[k]) then begin
	    {NOTE 02jan04: dont know what the "nre" was supposed to be... added n_a[icp]=n_a[i] to the condition below}
	    bx[k]:=(eqatomrec(n_para[i,j],n_para[icp,k]) and (n_a[icp,k,1]=n_a[i,j,1]));
	    if bx[k] then begin
	       b[j]:=true;
	       if (n_a[icp,k,1]<>n_a[i,j,1]) then begin
		  nre:=nre+1;
	       end;
	       {write('DBG ',i:1,' ',j:1,' ',k:1,'   ');
	       for o:=1 to 3 do write(n_a[i,j,o]:5:1); write('  ');
	       for o:=1 to 3 do write(n_a[icp,k,o]:5:1); writeln;}
	       nm:=nm+1;
	    end;
	 end;
      end;
      for j:=1 to nan[i] do if not(b[j]) then begin {new flipped atoms}
	 {writeln(j:1,' ',n_a[i,j,1]:1);}
	 n_1:=n_1+1;
	 po_1[n_1,1]:=i; po_1[n_1,2]:=j;
	 for o:=1 to 3 do r_1[n_1,o]:=bound(n_r[i,j,o]+rt[i,o]);
	 t_1[n_1]:=n_a[i,j,1];
      end;
      for j:=1 to nan[icp] do if not(bx[j]) then begin {old flipped atoms}
	 {writeln(j:1,' ',n_a[k,j,1]:1);}
	 n_0:=n_0+1;
	 for o:=1 to 3 do
	    r_0[n_0,o]:=bound(n_r[icp,j,o]+rt[icp,o]); {note rt[i]==rt[icp]}
	 t_0[n_0]:=n_a[icp,j,1];
      end;
      for j:=1 to nan[i] do if b[j] then begin {identic atoms}
	 n_2:=n_2+1;
	 for o:=1 to 3 do
	    r_2[n_2,o]:=bound(n_r[i,j,o]+rt[i,o]);
	 t_2[n_2]:=n_a[i,j,1];
      end;
   end else begin {flipped node(s): dont compare!}
      i:=fln[i0];
      for j:=1 to nan[i] do begin
	 n_1:=n_1+1; {new atoms}
	 po_1[n_1,1]:=i; po_1[n_1,2]:=j;
	 for o:=1 to 3 do r_1[n_1,o]:=bound(n_r[i,j,o]+rt[i,o]);
	 t_1[n_1]:=n_a[i,j,1];
      end;
      icp:=tnodes+i0;
      for j:=1 to nan[icp] do begin
	 n_0:=n_0+1; {old atoms}
	 for o:=1 to 3 do r_0[n_0,o]:=bound(n_r[icp,j,o]+rt[icp,o]);
	 t_0[n_0]:=n_a[icp,j,1];
      end;
      if(debug[3]) then writeln('cmpchl : ',nfl:1,' ',i0:1,' ',fln[i0]:1,' ',nan[i]:1,' ',nan[icp]:1,' (nfl i fln[i] nan[i] nan[icp])');
   end;
end; { CMPCHL }
procedure CMPCHLA(swo :integer; flst : flist; nfl : integer;fln:flist; Var nm:integer );
{
First eliminate atoms that stayed unchanged upon flip.
Then associate each flipped atom with an old atom,
determine chemistry of a new atom. Algorithm:
create lists of atoms of the same kind (comparing
last 3 records of n_para). Then randomly pair old-new
atoms of the "same kind" (ei the same orbit, controlled
by the same parameters).
}
Var
   i0,j0,i,j,k,o,ifl : integer;
   nre,icp,nl,nl0    : integer;
   por		     : myreal;
   b,bx		     : array[1..npth,1..at_of_node] of boolean;
   bfl		     : array[1..npth] of boolean;
   l		     : array[1..npth] of i2;
   procedure chassign(n : integer);
   Var j0,k,icp,o : integer;
   begin
      {note: i0, i, j are protected!}
      j0:=l[n,1]; k:=l[n,2];
      icp:=tnodes+j0;
      b[i0,j]:=true; bx[j0,k]:=true;
      n_a[i,j,1]:=n_a[icp,k,1];
      para2xyz(n_para[i,j],n_a[i,j,1],1,1,n_r[i,j]); {new position is chem-dep.}
      {old atoms}
      n_0:=n_0+1;
      for o:=1 to 3 do r_0[n_0,o]:=bound(n_r[icp,k,o]+rt[icp,o]);
      t_0[n_0]:=n_a[icp,k,1];
      {new atoms}
      n_1:=n_1+1;
      po_1[n_1,1]:=i; po_1[n_1,2]:=j;
      for o:=1 to 3 do r_1[n_1,o]:=bound(n_r[i,j,o]+rt[i,o]);
      t_1[n_1]:=n_a[i,j,1];
   end; { chassign }
begin
   n_0:=0; n_1:=0; n_2:=0; nm:=0;
   for i:=1 to nfl do begin
      ifl:=fln[i];
      {the "else" condition takes into account case couple=1}
      {if(swo=1) then bfl[i]:=false else begin}
	 bfl[i]:=false;
	 for j:=1 to flst[0] do if (fln[i]=flst[j]) then bfl[i]:=true;
	 {bfl[i]:=((fln[i]=flst[1]) or (fln[i]=flst[2]));}
      {end;}
      for j:=1 to nan[ifl] do b[i,j]:=false;
      for j:=1 to nan[i+tnodes] do bx[i,j]:=false;
   end;
   nre:=0;
   {Stage 1 (eliminate identic atoms):
   compare node by node, ommit flipped nodes}
   for i0:=1 to nfl do if not(bfl[i0]) then begin
      i:=fln[i0];
      icp:=tnodes+i0;
      for j:=1 to nan[i] do begin {loop new atoms}
	 for k:=1 to nan[icp] do if not(bx[i0,k]) then begin {old atoms}
	    bx[i0,k]:=eqatomrec(n_para[i,j],n_para[icp,k]);
	    if bx[i0,k] then begin
	       b[i0,j]:=true;
	       {in combination with LG, chemistry may have been whatever}
	       n_a[i,j,1]:=n_a[icp,k,1];
	       {displacements are chemically dependent, so need copy too}
	       n_r[i,j]:=n_r[icp,k]; 
	    end;
	 end;
      end;
      for j:=1 to nan[i] do if b[i0,j] then begin {identic atoms}
	 n_2:=n_2+1;
	 for o:=1 to 3 do
	    r_2[n_2,o]:=bound(n_r[i,j,o]+rt[i,o]);
	 t_2[n_2]:=n_a[i,j,1];
      end;
   end;
   {Stage 2: take an atom; scan the list of old atoms to create
   a sublist of alike old atoms; randomly choose one of them, and
   set chemistry of the new atom accordingly}
   for i0:=1 to nfl do begin {list of new atoms}
      i:=fln[i0];
      for j:=1 to nan[i] do if not(b[i0,j]) then begin {loop new atoms}
	 nl:=0;
	 for j0:=1 to nfl do begin
	    icp:=tnodes+j0;
	    for k:=1 to nan[icp] do if not(bx[j0,k]) then begin {old atoms}
	       if (n_a[i,j,2]=n_a[icp,k,2]) then begin
		  nl:=nl+1;
		  l[nl,1]:=j0; l[nl,2]:=k;
	       end;
	    end;
	 end;
	 nl0:=nl;
	 if (nl=1) then chassign(nl);
	 if (nl>1) then begin
	    por:=randu(seed)*nl; nl0:=round(por+0.5);
	    chassign(nl0);
	 end;
	 {writeln(i0:1,' ',j:1,' ',n_a[i,j,2]:1,' ',nl:1,' -> ',nl0:1);}
      end;
   end;
   {Stage 3: randomly reassign chemistries}
   for i0:=1 to nfl do begin {list of new atoms}
      i:=fln[i0];
      for j:=1 to nan[i] do if not(b[i0,j]) then begin {loop new atoms}
	 nl:=0;
	 for j0:=1 to nfl do begin
	    icp:=tnodes+j0;
	    for k:=1 to nan[icp] do if not(bx[j0,k]) then begin {old atoms}
	       nl:=nl+1;
	       l[nl,1]:=j0; l[nl,2]:=k;
	    end;
	 end;
	 nl0:=nl;
	 if (nl=1) then chassign(nl);
	 if (nl>1) then begin
	    por:=randu(seed)*nl; nl0:=round(por+0.5);
	    chassign(nl0);
	 end;
	 {writeln(i0:1,' ',j:1,' ',n_a[i,j,2]:1,' ',nl:1,' -> ',nl0:1);}
      end;
   end;
   nm:=0;
   for i:=1 to nfl do
      for j:=1 to nan[fln[i]] do if not(b[i,j]) then nm:=nm+1;
   if ((nm>0) and not(chvar)) then begin
      writeln('Couldnt register atoms before/after flip for node ',flst[1]:1,' :');
      writeln(nflp:1,' ',ntflp:1,' cmpchla: ',nm:1,' atoms left out.');
      writeln('Exporting copa.err.');
      export_config_fl('copa.err');
      halt;
   end;
end;{CMPCHLA}

procedure ieinit; {loop over images of the unit cell}
Var
   p1,p2,p3 : integer;
   ar	    : myreal;
   w,y	    : v3;
   ide	    : boolean;
begin
   npx:=0;
   for p1:=-iext[1] to iext[1] do
      for p2:=-iext[2] to iext[2] do
	 for p3:=-iext[3] to iext[3] do begin
	    ide:=((p1=0) and (p2=0) and (p3=0));
	    y[1]:=p1; y[2]:=p2; y[3]:=p3;
	    ortho3(y,w);
	    ar:=sqrt(w[1]*w[1]+w[2]*w[2]+w[3]*w[3]);
	    if ((ar<rcut) and not(ide)) then begin
	       npx:=npx+1;
	       if(npx>iucpp)then begin
		  writeln('utabinit (tile_mc): increase const iucpp.');
		  halt;
	       end;
	       px[npx]:=w; {!! in [A] !!}
	    end;
	 end;
end; { ieinit }
function add_if(x : integer; ql:integer;q:vpath; Var v : flist; Var nv:integer):integer;
Var
   i   : integer;
   add : boolean;
begin
   add:=true;
   add_if:=0;
   for i:=1 to ql do if(x=q[i]) then add:=false;
   if(add) then for i:=1 to nv do if(x=v[i]) then add:=false;
   if(add) then begin
      nv:=nv+1;
      v[nv]:=x;
      add_if:=1;
   end;
end; { add_if }
procedure boundary(iss,dp,sg: integer; p,q:vpath; var v:flist; var nv:integer);
var
   ref,ref1,j0	   : integer;
   i,j,k,k0,k1,l,n : integer;
   nv0,nv1,nn	   : integer;
   b		   : array[1..t_flip] of boolean;
begin
   {writeln(' --DBG-BOUNDARY');}
   nv:=0;
   i:=0; j:=iss; ref:=p[1]; p[1]:=0; repeat i:=i+1;
      if((i>1) and (i<dp)) then begin
	 k0:=po(ref+5); {previous vertex dir}
	 k1:=po(ref+p[i]); {following vertex dir}
	 {writeln(' next ',iss:1,' ',j:1,' ',k0:1,' ',k1:1);}
	 k:=k1; repeat
	    k:=k+sg;
	    l:=lt[j,po(k)];
	    {write(po(k):1,':',l:1,' ');}
	    if(l>0) then n:=add_if(l,dp,q,v,nv);
	 until (po(k+sg)=k0);
	 {writeln;}
      end;
      ref:=po(ref+p[i]);
      j:=lt[j,ref];
   until (i=dp);
   nv0:=1; nv1:=nv; repeat nn:=0;
      for i:=nv0 to nv1 do begin
	 for j:=1 to 10 do if lt[v[i],j]>0 then nn:=nn+add_if(lt[v[i],j],dp,q,v,nv);
      end;
      if(nn>0) then begin
	 nv0:=nv1+1;
	 nv1:=nv+nn;
      end;
   until (nn=0);
end; { boundary }
function findop(a1,a2 : i2):integer;
Var
   found : boolean;
   i,j,o : integer;
begin
   found:=false; findop:=0;
   for i:=1 to n_sym do if not(found) then begin
      if((sym[i,a1[1]]=a2[1]) and (sym[i,a1[2]]=a2[2])) then begin
	 findop:=i; found:=true;
      end;
   end;
end; { findop }
procedure ltupif(n1,n2,x : integer; Var ne:integer );
Var
   j1,j2,i,j,o,k,y2,p,pk,q : integer;
   found		   : boolean;
begin
   found:=false;
   for j1:=-1 to 1 do if not(found) then
      for j2:=-1 to 1 do begin if not(found) then
	 y2:=0;
	 for o:=1 to 5 do begin if not(found) then
	    k:=bt[n1,o]-bt[n2,o]-j1*bai[1,o]-j2*bai[2,o];
	    if (k*k=1) then begin p:=o; pk:=k end;
	    y2:=y2+k*k;
	 end;
	 if((y2=1) and not(found)) then begin
	    found:=true;
	    q:=(p-1)*2+1;
	    if(pk=-1) then q:=po(q+5);
	    ne:=ne+1;
	    lt[n2,q]:=n1*x;
	    q:=po(q+5);
	    lt[n1,q]:=n2*x;
	    if debug[2] then begin
	       if(x=1) then writeln('  created linkage ',n1:1,'-',n2:1);
	       if(x=0) then writeln('  erased linkage ',n1:1,'-',n2:1);
	    end;
	 end;
      end;
end;
procedure ltupif3(n1,n2,x : integer; Var ne:integer);
Var
   j1,j2,i,j,o,k,y2,p,pk,q : integer;
   k1,k2,kk		   : integer;
   found		   : boolean;
begin
   k1:=((n1-1) div tnodes0);
   k2:=((n2-1) div tnodes0);
   kk:=k1-k2;
   if(kk>1) then kk:=kk-multi;
   if(kk<-1) then kk:=kk+multi;
   if(kk*kk<2) then begin
      found:=false;
      for j1:=-1 to 1 do if not(found) then
	 for j2:=-1 to 1 do begin if not(found) then
	    y2:=0;
	    for o:=1 to 5 do begin if not(found) then
	       k:=bt[n1,o]-bt[n2,o]-j1*bai[1,o]-j2*bai[2,o];
	       if (k*k=1) then begin p:=o; pk:=k end;
	       y2:=y2+k*k;
	    end;
	    if((y2=1) and not(found)) then begin
	       found:=true;
	       ne:=ne+1;
	       q:=(p-1)*2+1;
	       if(pk=-1) then q:=po(q+5);
	       ltm[n2,kk,q]:=n1*x;
	       q:=po(q+5);
	       ltm[n1,-kk,q]:=n2*x;
	       {writeln('y2=1 ',n1:1,' ',n2:1,' ',kk:1);}
	       if (debug[2] and (kk<>0)) then begin
		  if(x=1) then writeln('  created inter-lay-linkage ',n1:1,'-',n2:1,' height=',kk:1);
		  if(x=0) then writeln('  erased inter-lay-linkage ',n1:1,'-',n2:1,' height=',kk:1);
	       end;
	    end;
	    if((y2=0) and not(found) and (kk<>0)) then begin
	       found:=true;
	       {writeln('y2=0 ',n1:1,' ',n2:1,' ',kk:1);}
	       ltm[n2,kk,0]:=n1*x;
	       ltm[n1,-kk,0]:=n2*x;
	       if debug[2] then begin
		  if(x=1) then writeln('  created z-link ',n1:1,'-',n2:1);
		  if(x=0) then writeln('  erased z-link ',n1:1,'-',n2:1);
	       end;
	    end;
	 end;
   end;
end;
function flippable(i  : integer):boolean; begin flippable:=(nv[i]=3) end;
procedure check_flp(swo,i,cpl:integer; Var x:flist; Var flpo,nn_mu:integer; Var nn_list:flist);
{Always flip hexagon (3-fold vertex); in case a fat hexagon (nu=2)
could be flipped if a neighboring skinny hex. flipped, perform
DOUBLE flip, skinny first and then the fat. The incomplete
fat hexagon is recognized from 4-fold vertex with nu=2, it has
to have 2 neighbors with nu=2 and two with nu=1. Moreover,
one of the two nu=1 neighboring vertices has to be 3-fold
vertex (flippable skinny hex.)
Variable x caries information: x[0]=0 if the vertex
is not flippable; x[0]=1 and x[1]=i if it's 3-fold vertex;
and x[0]=2, x[1]=jf and x[2]=i for "double" flip.}
{NEW!!!
first parameter swo=1 means flip_object!!!
}
{NEW july 2005 - generalized flip}
Var
   nuii,o,p,q,k,j,l1,l2,pl,io : integer;
   n1,n2,jf,sg,pse,fl0,fl1,ii : integer;
   done,flable,fpab	      : boolean;
   xr			      : myreal;
   ap			      : vpath;
   nli,nvb,ix,m1,nvb3	      : integer;
   li			      : flist;
   x5,y5		      : v5;
   ix5			      : i5;
   x1,x2		      : v10;
   w			      : v3;
   dd			      : myreal;
   zflippable		      : boolean;
   ap3			      : array[-1..1] of vpath;
   li3			      : array[-1..1] of flist;
   nedg			      : i2;
   ipthpo		      : integer;
begin
   {dump objects}
   {for j:=1 to tnodes do
      for k:=1 to n_objmu[j] do begin
	 l1:=n_objpo[j,k,1];
	 writeln(j:1,' ',k:1,' ',obj_name[l1]);
      end;
   halt;}
   for j:=0 to t_flip do x[j]:=0;
   nn_mu:=0; for j:=1 to t_flip do nn_list[j]:=0;
   flpo:=0;
   if((swo=1) and (n_objmu[i]>0)) then begin {case flip-object mode}
      xr:=randu(seed)*n_objmu[i];
      pl:=round(xr+0.5);
      if(pl=0)then pl:=1; if(pl>n_objmu[i])then pl:=n_objmu[i];
      l1:=n_objpo[i,pl,1]; {object}
      l2:=n_objpo[i,pl,2]; {object_orientation}
      if(flip_man) then begin
	 k:=0;
	 for j:=1 to n_objmu[i] do begin
	    l1:=n_objpo[i,j,1]; {object}
	    l2:=n_objpo[i,j,2]; {object_orientation}
	    if((obj_nfl[l1]>1) or (nobcolors[l1]>1)) then begin k:=k+1; x[k]:=j; end;
	 end;
	 if(k>0) then begin
	    {READING STDINP DOES NOT WORK - ITS ALL ZEROES...}
	    if(k>1) then begin
	       write('oflip-manual: there are ',k:1,' flippable objects - choose one: '); readln(o);
	    end else begin
	       o:=1;
	    end;
	    l1:=n_objpo[i,x[o],1]; {object}
	    l2:=n_objpo[i,x[o],2]; {object_orientation}
	    writeln('... flipping object ',o:1);
	 end else begin
	    writeln('Node ',i:1,' has no flippable objects.');
	    halt;
	 end;
	 x[0]:=0; x[1]:=0;
      end;
      {writeln('xxxdbg ',l1:1,' ',nobcolors[l1]:1);}
      {writeln('yyydbg ',obj_nfl[l1]:1);}
      if(obj_nfl[l1]*nobcolors[l1]>1) then begin
	 flpo:=l1;
	 if(obj_nipth[l1]>0) then begin
	    xr:=randu(seed)*obj_nipth[l1];
	    o:=round(xr+0.5);
	    ipthpo:=obj_ipth[l1,o];
	 end else begin
	    ipthpo:=0;
	 end;
	 {writeln('DBGff object ',obj_name[l1],' ',obj_nipth[l1]:1,' choice=',o:1,' ',ipthpo:1);}
	 {the old flip-object mode, when the object has atoms bound to it}
	 if(obj_nfl[l1]=1) then begin
	    x[0]:=1; x[1]:=i;
	    nn_mu:=1; nn_list[1]:=i;
	    copy_nrec(nn_list[1],tnodes+1);
	    copy_trec(nn_list[1],tnodes+1);
	    if deco then COPY_AREC(nn_list[1],tnodes+1);
	 end else begin {otherwise we try to reshuffle tiling}
	    if(debug[2]) then writeln('attempting general-flip ',obj_name[l1],' from root-vert ',i:1);
	    nvb:=obj_pth[l1];
	    {pse determines sense in which the path closes}
	    pse:=0; for j:=2 to nvb do pse:=pse+obj_pth_def[l1,l2,j];
	    if(pse>0) then pse:=1 else pse:=-1;
	    fpab:=fpath(i,nvb,obj_pth_def[l1,l2],ap);
	    for j:=1 to nvb do begin
	       if (ap[j]=0) then begin
		  writeln('d04 General flip: inconsistent object def - probably need to increase ulevel');
		  halt;
	       end;
	    end;
	    boundary(i,nvb,pse,obj_pth_def[l1,l2],ap,li,nli);
	    {writeln('DBG after bound nli=',nli:1);}
	    zflippable:=true;
	    if(io3d>0) then begin
	       for o:=-1 to 1 do if (o<>0) then begin
		  for j:=1 to nvb do if(ltm[ap[j],o,0]=0) then begin
		     zflippable:=false;
		  end else begin
		     ap3[o,j]:=ltm[ap[j],o,0];
		  end;
	       end;
	    end;
	    if(zflippable) then begin
	       if(io3d>0) then begin
		  for o:=-1 to 1 do if (o<>0) then begin
		     ix:=ltm[i,o,0];
		     boundary(ix,nvb,pse,obj_pth_def[l1,l2],ap3[o],li3[o],nli);
		  end;
	       end;
	       for j:=1 to nli do begin {add vertices inside boundary to node-update list}
		  x[0]:=x[0]+1; x[x[0]]:=li[j];
		  nn_mu:=nn_mu+1; nn_list[nn_mu]:=li[j];
	       end;
	       for j:=1 to nvb do begin {add boundary to node-update list}
		  nn_mu:=nn_mu+1;
		  nn_list[nn_mu]:=ap[j];
	       end;
	       m1:=-1; if multi=2 then m1:=1;
	       if(io3d>0) then for o:=m1 to 1 do if(o<>0) then begin
		  for j:=1 to nli do begin {add vertices inside boundary to node-update list}
		     nn_mu:=nn_mu+1; nn_list[nn_mu]:=li3[o,j];
		  end;
		  for j:=1 to nvb do begin {add boundary to node-update list}
		     nn_mu:=nn_mu+1;
		     nn_list[nn_mu]:=ap3[o,j];
		  end;
	       end;
	       if debug[2] then begin
		  writeln(' check_flp ',obj_name[l1],' root-vert ',i:1,' ');
		  write('  boundary nvb=',nvb:1,' : ');
		  for j:=1 to nvb do begin
		     write(ap[j]:1,' ');
		  end; writeln;
		  write('  internal-verts nli=',nli:1,' : ');
		  for j:=1 to nli do begin
		     write(li[j]:1,' ');
		  end; writeln;
	       end;
	       fl0:=n_obj_fl[i,pl];
	       repeat
		  xr:=randu(seed)*obj_nfl[l1]; n1:=round(xr+0.5);
		  if(n1=0) then n1:=1; if(n1>obj_nfl[l1]) then n1:=obj_nfl[l1];
	       until (n1<>fl0);
	       fl1:=n1;
	       {write(obj_nfl[l1]:1,'  ',fl0:1,' ',n1:1,'   ');
	       write(obj_flb[l1,l2,fl0,1]:1,':',obj_flb[l1,l2,fl0,2]:1,'   ');
	       writeln(obj_flb[l1,l2,fl1,1]:1,':',obj_flb[l1,l2,fl1,2]:1,'   ');}
	       io:=findop(obj_flb[l1,l2,fl0],obj_flb[l1,l2,fl1]);
	       if(io=0) then begin
		  writeln('generalized flip error: invalid symmetry operation');
		  write('  fl0->');
		  for o:=1 to 2 do write(obj_flb[l1,l2,fl0,o]:1,' ');
		  write('  fl1->');
		  for o:=1 to 2 do write(obj_flb[l1,l2,fl1,o]:1,' ');
		  writeln;
		  halt;
	       end;
	       if debug[2] then writeln('  performing symop # ',io:1,'  ',fl0:1,'->',fl1:1);
	       {save current state}
	       for j:=1 to nn_mu do copy_nrec(nn_list[j],tnodes+j);
	       for j:=1 to nn_mu do copy_trec(nn_list[j],tnodes+j);
	       if deco then for j:=1 to nn_mu do
		  COPY_AREC(nn_list[j],tnodes+j);
	       {for pat-instances, we need to clear all old ones BEFORE update move}
	       if(nipatot>0) then for j:=1 to nn_mu do new_pat_instance(nn_list[j],1);
	       {update connectivity}
	       for o:=1 to 2 do nedg[o]:=0;
	       if(nli>0) then begin
		  for j:=1 to nvb do begin
		     for k:=1 to nli do begin
			n1:=ap[j];
			n2:=li[k];
			ltupif(n1,n2,0,nedg[1]);
		     end;
		  end;
		  {first nli entries in nn_list are vertices at cur level inside boundary. Next nvb is boundary at this level, rest are up/down levels.}
		  if(io3d>0) then begin
		     for j:=nli+1 to nn_mu do begin
			for k:=1 to nli do begin
			   n1:=nn_list[j];
			   n2:=li[k];
			   ltupif3(n1,n2,0,nedg[1]);
			end;
		     end;
		  end;
		  {update vertices inside.boundary at this level in between themselves}
		  for j:=1 to nli-1 do begin
		     for k:=j+1 to nli do begin
			n1:=li[j];
			n2:=li[k];
			ltupif(n1,n2,0,nedg[1]);
		     end;
		  end;
	       end;
	       for j:=1 to nli do begin
		  ii:=li[j];
		  for o:=1 to 5 do x5[o]:=bt[ii,o]-obj0[l1,l2,o]-bt[i,o];
		  bound5r(x5,w);
		  for o:=1 to 5 do y5[o]:=0.0;
		  for o:=1 to 5 do begin
		     n1:=sym[io,(o-1)*2+1];
		     for p:=1 to 5 do y5[p]:=y5[p]+bs[n1,p]*x5[o];
		  end;
		  if (debug[2]) then begin
		     writeln('vertex ',ii:1,' : ');
		     for o:=1 to 5 do write(x5[o]:1:2,' '); writeln(' old');
		     for o:=1 to 5 do write(y5[o]:1:2,' '); writeln(' new');
		  end;
		  for o:=1 to 5 do ix5[o]:=round(y5[o]+bt[i,o]+obj0[l1,l2,o]);
		  bound5(ix5,w);
		  bt[ii]:=ix5; {new vertex position}
		  nu[ii]:=nui(ix5);
		  rt[ii]:=w;
		  rt[ii,3]:=((ii-1) div tnodes0)/multi;
		  if(ipthpo<>0) then begin
		     tvl[ii]:=ipthpo;
		     {writeln('DBG vertex ',ii:1,' newtvl=',ipthpo:1);}
		  end;
	       end;
	       for k:=1 to nli do begin
		  for o:=1 to 10 do lt[li[k],o]:=0;
	       end;
	       for j:=1 to nvb do begin
		  for k:=1 to nli do begin
		     n1:=ap[j];
		     n2:=li[k];
		     ltupif(n1,n2,1,nedg[2]);
		  end;
	       end;
	       if(nli>0) then begin
		  if(io3d>0) then begin
		     for j:=nli+1 to nn_mu do begin
			for k:=1 to nli do begin
			   n1:=nn_list[j];
			   n2:=li[k];
			   ltupif3(n1,n2,1,nedg[2]);
			end;
		     end;
		  end;
		  for j:=1 to nli-1 do begin
		     for k:=j+1 to nli do begin
			n1:=li[j];
			n2:=li[k];
			ltupif(n1,n2,1,nedg[2]);
		     end;
		  end;
	       end;
	       for j:=1 to tnodes do fl[j]:=0; fl[i]:=1;
	       if(nedg[1]<>nedg[2]) then begin
		  {writeln('error: number of tiling edges not conserved upon general flip');
		  export_config_fl('copa.tmp');
		  halt;}
	       end;
	       {halt;}
	    end else begin
	       x[0]:=0;
	    end;
	 end;
      end else begin
	 x[0]:=0; x[1]:=0;
      end;
   end else if (swo=0) then begin {case flip-vertex mode}
      if (nv[i]=3) then begin x[0]:=1; x[1]:=i end else if (cpl=1) then begin
	 if (nv[i]=4) then begin
	    nuii:=dperp_offset(-nu_offs,nu[i]);
	    if (abs(nuii)=2) then begin
	       j:=0; n2:=0; n1:=0; jf:=0; done:=false; repeat j:=j+1;
		  if (lt[i,j]>0) then begin
		     k:=dperp_offset(-nu_offs,nu[lt[i,j]]);
		     if (abs(k)=2) then n2:=n2+1;
		     if (abs(k)=1) then n1:=n1+1;
		     if ((abs(k)=1) and (nv[lt[i,j]]=3)) then jf:=j;
		  end;
	       until (j=10);
	       if ((n1=2) and (n2=2) and (jf>0)) then begin
		  x[0]:=2; x[1]:=lt[i,jf]; x[2]:=i;
	       end;
	    end;
	 end;
      end;
   end else if (swo=2) then begin {instance flip}
      flable:=false;
      for k:=1 to n_objmu[i] do begin
	 l1:=n_objpo[i,k,1]; {object}
	 l2:=n_objpo[i,k,2]; {object_orientation}
	 if(obj_nipat[l1]>0) then flable:=true;
	 {writeln('xxx ',l1:1,' ',l2:1,' ',obj_nipat[l1]:1);}
      end;
      if (flable) then begin
	 nn_mu:=1; nn_list[1]:=i;
	 for k:=1 to 10 do if (lt[i,k]>0) then begin
	    nn_mu:=nn_mu+1;
	    nn_list[nn_mu]:=lt[i,k];
	 end;
	 {save}
	 for j:=1 to nn_mu do copy_nrec(nn_list[j],tnodes+j);
	 for j:=1 to nn_mu do copy_trec(nn_list[j],tnodes+j);
	 if deco then for j:=1 to nn_mu do
	    COPY_AREC(nn_list[j],tnodes+j);
	 x[0]:=1; x[1]:=i;
	 if(instances) then begin
	    for j:=1 to nn_mu do new_pat_instance(nn_list[j],0);
	    new_pat_instance(i,0);
	 end;
      end else begin
	 x[0]:=0; x[1]:=0;
      end;
   end;
   {we did not need to update nv, we do it once here}
   for j:=1 to tnodes do begin
      k:=0; for o:=1 to 10 do if(lt[j,o]>0) then k:=k+1;
      nv[j]:=k;
   end;
end;
procedure check_sigmaflp(ii,ulv:integer; Var nnb:integer; Var blist:flist; Var flipped:boolean);
{
 FLIP "Sigma"-vertex IN A GENERAL TWO-LEVEL TILING
}
Const
   bmax	 = 100;
   nsurf = 4;
Var
   o,k,i,j,l,mydbg   : integer;
   nd,ed,ltp,d1r,d2r : integer;
   vini,vcur,vnew    : integer;
   dini,dcur,dnew    : integer;
   d,n		     : i2;
   x		     : i4;
   xio		     : array[-1..1] of i4;
   bclosed	     : boolean;
   fv		     : i5;
   w		     : v3;
   zflippable	     : boolean;
   lz		     : array[-1..1] of integer;
   x4s		     : array[1..4,-1..1] of i10;
   x2		     : array[-1..1] of integer;
   upd_fast	     : boolean;
   nnu		     : integer;
   ulist	     : flist;
   procedure Zstuff;
   Var o,l,i,j,k,p1 : integer;
   begin
      l:=nnb; {5+boundary; this tiling nodelist already set}
      for o:=-1 to 1 do lz[o]:=0; 
      for o:=-1 to 1 do if ((o<>0) and (zflippable)) then begin
	 if (node_z_alligned(ii,o)<>ltm[ii,o,0]) then begin
	    writeln('Zstuff LTM problem: ',ltm[ii,o,0]:1,' ',node_z_alligned(ii,o):1);
	    halt;
	 end;
	 lz[o]:=ltm[ii,o,0]; {now implemented directly into the tables}
	 if (lz[o]=0) then begin {the top/bottom non-z-alligned vertex}
	    if (ltc[ii,o]=0) then begin
	       writeln('Zstuff VC error ',ii:1,': ltm=0 AND ltc=0.');
	       halt;
	    end;
	    l:=l+1; blist[l]:=ltc[ii,o];
	    x2[o]:=ltc[ii,o];
	    p1:=lt[ltc[ii,o],ed];
	    if(p1>0) then begin
	       l:=l+1; blist[l]:=p1;
	    end;
	 end else begin
	    l:=l+1; blist[l]:=lz[o];
	    x2[o]:=lz[o];
	    p1:=lt[lz[o],po(ed+5)];
	    if(p1>0) then begin
	       l:=l+1; blist[l]:=p1;
	    end;
	 end;
	 for k:=1 to nsurf do if (zflippable) then begin
	    xio[o,k]:=ltm[x[k],o,0]; {pointer to node}
	    if(xio[o,k]>0) then begin
	       l:=l+1; blist[l]:=xio[o,k];
	    end else zflippable:=false;
	 end;
      end;
      if (flip_man and not(zflippable)) then begin
	 writeln('Vertex ',ii:1,' is not z-flippable.');
      end;
      nnb:=l;
   end; { Zstuff }
   procedure Zupdate;
   Var
      o,i,j,k  : integer;
      p1,p2,p3 : integer;
   begin
      {write('BEF ',lz[-1]:1,' ',lz[1]:1);}
      for o:=-1 to 1 do if (o<>0) then begin
	 for j:=1 to 2 do begin {erase old bonds}
	    p1:=xio[o,j]; p2:=d[j];
	    if (ltm[p1,-o,po(p2+5)]<>ii) then begin
	       writeln('Zupdate: broken VC: ',ii:1,'<>',ltm[p1,-o,po(p2+5)]:1,', p1=',p1:1,'(->tiling.err)');
	       export_tiling('tiling.err');
	       halt;
	    end;
	    ltm[ii,o,p2]:=0;
	    nvm[ii,o]:=nvm[ii,o]-1;
	    ltm[p1,-o,po(p2+5)]:=0;
	    nvm[p1,-o]:=nvm[p1,-o]-1;
	 end;
	 for k:=1 to 2 do begin {create new bonds}
	    l:=xio[o,k+2];
	    if(k=1) then j:=2 else j:=1;
	    p2:=d[j];
	    ltm[l,-o,p2]:=ii;
	    nvm[l,-o]:=nvm[l,-o]+1;
	    ltm[ii,o,po(p2+5)]:=l;
	    nvm[ii,o]:=nvm[ii,o]+1;
	 end;
	 {write('CHECK '); for j:=1 to 10 do write(ltm[ii,o,j]:1,' '); writeln;}
	 p1:=ltm[ii,o,0];
	 if (p1>0) then begin {ii and p1 were alligned: erase}
	    ltm[ii,o,0]:=0;
	    ltm[p1,-o,0]:=0;
	    ltc[ii,o]:=p1;
	    ltc[p1,-o]:=ii;
	    {NEW STUFF: the vertex above the flipping vertex is 3-fold}
	    p2:=ltm[ii,o,po(ed+5)];
	    if(p2>0) then begin
	       p3:=ltm[p2,-o,ed];
	       if(p3<>ii) then begin
		  writeln('Zupdate: error in ltm tables: ',ii:1,'!=',p3:1,'(p1=',p1:1,', p2=',p2:1,'edi=',po(ed+5):1,')');
		  halt;
	       end;
	       ltm[ii,o,po(ed+5)]:=0;
	       ltm[p3,-o,ed]:=0;
	    end;
	 end else begin {ii and p1 were NOT alligned: create z-allignement entry}
	    {p1:=node_z_alligned(ii,o);}
	    p1:=ltc[ii,o];
	    ltm[ii,o,0]:=p1;
	    ltm[p1,-o,0]:=ii;
	    ltc[ii,o]:=0;
	    ltc[p1,-o]:=0;
	    {NEW STUFF: the vertex above the flipping vertex is 3-fold}
	    p2:=lt[p1,ed];
	    if(p2>0) then begin
	       ltm[ii,o,ed]:=p2;
	       ltm[p2,-o,po(ed+5)]:=ii;
	    end;
	 end;
      end;
      for o:=-1 to 1 do if (o<>0) then begin
	 lz[o]:=ltm[ii,o,0]; {update node z-allignments}
      end;
      {writeln('   AFT ',lz[-1]:1,' ',lz[1]:1);}
   end; { Zupdate }
begin
   upd_fast:=false;
   flipped:=false;
   mydbg:=0;
   for o:=1 to 4 do x[o]:=0; d[1]:=0; d[2]:=0;
   if(nv[ii]=2) then begin
      nd:=0;
      for o:=1 to 10 do if (lt[ii,o]>0) then begin
	 nd:=nd+1; d[nd]:=o; x[nd]:=lt[ii,o];
      end;
      ed:=0;
      if(d[2]-d[1]=6) then begin
	 ed:=po(d[1]+3);
      end else if(d[2]-d[1]=4) then begin
	 ed:=po(d[1]+2+5);
      end else begin
	 {writeln("check_sigmaflp : bad 2-fold vertex (doesnt divide angle as 6+4)");
	 halt;}
      end;
      if(ed>0) then begin
	 if((lt[x[1],ed]>0) and (lt[x[2],ed]>0)) then begin
	    x[3]:=lt[x[1],ed]; x[4]:=lt[x[2],ed];
	 end;
      end;
   end;
   if((x[1]>0) and (x[2]>0) and (x[3]>0) and (x[4]>0)) then begin {FLIPOK}
      zflippable:=true;
      if ((io3d>0) and (multi>1)) then Zstuff; {returns zflippable}
      if(zflippable) then begin
      {boundary: take maximal turns to the right, and stop after reaching vertex1}
      {we start checking from direction opposite to the one along we moved last time}
      {note: sense of path bending doesnt matter}
	 vini:=x[1]; dini:=po(d[1]+5);
	 vcur:=vini; dcur:=dini;
	 nnb:=2; blist[1]:=ii; blist[2]:=vcur;
	 j:=0; repeat j:=j+1;
	    vnew:=0; dnew:=0;
	    for o:=dcur+1 to dcur+9 do begin
	       ltp:=lt[vcur,po(o)];
	       if(ltp>0) then if (ltp <> ii) then begin
	       vnew:=ltp; dnew:=po(o+5);
	       end;
	    end;
	    if(vnew>0) then begin vcur:=vnew; dcur:=dnew; end else begin
	       writeln('sigmaflip: could not find boundary contd ',vcur:1);
	       halt;
	    end;
	    bclosed:=(vcur=vini);
	    if not(bclosed) then begin nnb:=nnb+1; blist[nnb]:=vcur; end;
	 until (bclosed) or (j=bmax);
	 if((j=bmax) and not(bclosed)) then begin
	    writeln('flip boundary too long');
	    halt;
	 end;
	 {perform the flip}
	 fv:=bt[ii];
	 for o:=1 to 5 do fv[o]:=fv[o]+bs[d[1],o]+bs[ed,o]+bs[d[2],o];
	 bound5(fv,w);
	 {update connectivity tables}
	 if(upd_fast) then begin
	    {save}
	    for j:=1 to nnb do copy_nrec(blist[j],tnodes+j);
	    if deco then for j:=1 to nnb do COPY_AREC(blist[j],tnodes+j);
	    for j:=1 to nnb do copy_trec(blist[j],tnodes+j);
	    rt[ii]:=w;
	    bt[ii]:=fv; {new vertex position}
	    nu[ii]:=nui(fv);
	    d1r:=po(d[1]+5); d2r:=po(d[2]+5);
	    lt[ii,d[1]]:=0; lt[ii,d[2]]:=0;
	    lt[x[1],d1r]:=0; lt[x[2],d2r]:=0;
	    nv[x[1]]:=nv[x[1]]-1; nv[x[2]]:=nv[x[2]]-1;
	    lt[ii,d2r]:=x[3]; lt[ii,d1r]:=x[4];
	    lt[x[3],d[2]]:=ii; lt[x[4],d[1]]:=ii;
	    nv[x[3]]:=nv[x[3]]+1; nv[x[4]]:=nv[x[4]]+1;
	    if (io3d>0) then Zupdate;
	 end else begin
	    nnu:=0;
	    for i:=1 to 4 do begin
	       nnu:=nnu+1; ulist[nnu]:=x[i];
	       if((io3d>0) and (multi>1)) then for o:=-1 to 1 do if (o<>0) then begin
		  nnu:=nnu+1;
		  ulist[nnu]:=xio[o,i];
		  k:=nnu;
		  for j:=1 to 10 do if (lt[ulist[k],j]>0) then begin
		     nnu:=nnu+1; ulist[nnu]:=lt[ulist[k],j];
		  end;
		  if(ltc[ii,o]>0) then begin
		     nnu:=nnu+1; ulist[nnu]:=ltc[ii,o];
		     k:=nnu;
		     for j:=1 to 10 do if (lt[ulist[k],j]>0) then begin
			nnu:=nnu+1; ulist[nnu]:=lt[ulist[k],j];
		     end;
		  end;
	       end;
	    end;
	    nodes_union(nnu,ulist);
	    for i:=1 to nnu do begin
	       nnb:=nnb+1; blist[nnb]:=ulist[i];
	    end;
	    nodes_union(nnb,blist);
	    {save}
	    add_boundary(ulv-1,nnb,blist);
	    nodes_union(nnb,blist);
	    for j:=1 to nnb do copy_nrec(blist[j],tnodes+j);
	    rt[ii]:=w;
	    bt[ii]:=fv; {new vertex position}
	    nu[ii]:=nui(fv);
	    {writeln('copytrec ',nnb:1);
	    for i:=1 to nnb do write(blist[i]:1,' '); writeln;}
	    for j:=1 to nnb do copy_trec(blist[j],tnodes+j);
	    if deco then for j:=1 to nnb do COPY_AREC(blist[j],tnodes+j);
	    table_update(ii,nnu,ulist);
	 end;
	 flipped:=true;
      end else begin
	 for o:=1 to 4 do x[o]:=0; {NOT FLIPPABLE}
      end;  {zflippable}
   end;
end; { check_sigmaflp }
procedure mc_tiling(swo,cpl,vpin,lvl0,tlvl,kk,mcs,lab: integer; Var tpot,rate:myreal);
{
 swo : 2==[flip tiling node] 1==[flip s-broken object]
 cpl : flip coupling fat-skinny hexagon
 lvl0 : update level - how far to check what changed in the tiling upon flip
 kk : temperature cycle
 mcs: MC steps
 lab: label of "node_flip" or "object_flip"
 tpot: energy
 rate: acceptance rate
 }
Var
   i,j,k,l,o,upo,nuii,na0,na1,lvl   : integer;
   bf,fupd,fd			    : boolean;
   echnew,echold,por,acer,dchp	    : myreal;
   attempts,ii,jj,nla,nn0,nat,swptp : integer;
   u0,u1,u1x,u0x,dm,dux		    : myreal;
   du0,du2,dumt,x,y		    : myreal;
   force_reject,force_accept	    : boolean;
   try_reject			    : boolean;
   nn_mu,nn_mu0,n_surf,nn_mu1	    : integer;
   nn_tmu,nn_imu		    : integer;
   nn_list,nn_tlist,nn_ilist	    : flist;
   reject,flipped,swapok	    : boolean;
   flst				    : flist;
   bflst			    : array[1..t_flip] of boolean;
   dnz				    : array[0..z_ele] of integer;
   flipo			    : integer;
   nnma				    : integer;
   fff				    : nodelist;
begin
   export_err:=false;
   {writeln('export_err=',export_err);}
   {set swap type checking if orbits have mixed chemistry}
   swptp:=1; nla:=0;
   och_status; nat:=natoms;
   {NOTE: this test will incorrectly set swptp=2 for symmetry-broken orbits or objects, if those have mixed chemical occupancy}
   for i:=1 to maxorb do begin
      o:=0; for j:=1 to z_ele do if (naoch[i,j]>0) then begin
	 o:=o+1;
      end;
      if (o>1) then swptp:=2;
   end;
   if debug[3] then writeln('SWAP_TYPE : ',swptp:1);
   {check if there is a flippable object}
   lvl:=lvl0;
   {if(swo=1)then lvl:=1;}
   acer:=0.0;
   flipo:=0; {index of flipped object for general flip}
   fupd:=true;
   occupied_sites(natoc);
   for i:=1 to nnele do begin
      {writeln('DBG_CHVAR_BEG ',export_err,' ',i:1,' ',ele[i]:1,' ',mele_oc[i]:1);}
      {chnom[ele[i]]:=chcur[ele[i]];} {this must be set at the beginning only!}
   end;
   attempts:=mcs*tnodes+1;
   {attempts:=1;}
   if (flip_man or flip_replay) then begin
      dumt:=0.0;
      attempts:=nseq;
   end;
   for i:=1 to tnodes do fl[i]:=0;
   for jj:=1 to attempts do begin
      for k:=1 to tnodes do tvls[k]:=tvl[k];
      echold:=0.0; echnew:=0.0;
      if (flip_man or flip_replay) then begin
	 if (deco) then begin
	    {writeln('saving ',ivi:1);}
	    writeln(txta,jj:1,' ',vflip[jj]:1);
	    export_tiling_body(txta);
	    export_params(txta);
	    export_tiling_atoms(txta);
	 end;
	 ii:=vflip[jj];
	 if(sigmaflip) then begin
	    flst[0]:=0;
	    check_sigmaflp(ii,lvl,nn_mu,nn_list,flipped);
	    if(flipped) then begin flst[0]:=1; flst[1]:=ii; end;
	 end else begin
	    check_flp(swo,ii,cpl,flst,flipo,nn_mu,nn_list);
	    if (flst[0]=0) then begin
	       if (flip_replay or flip_man) then writeln('Error: vertex ',ii:1,' not flippable (nn=',nv[ii]:1,')');
	       halt;
	    end else begin
	    end;
	 end
      end else begin
	 repeat
	    eavnn[kk,0]:=eavnn[kk,0]+1.0;
	    eavnn[kk,lab]:=eavnn[kk,lab]+1.0;
	    {por:=randu(seed)*tnodes; ii:=round(por+0.5);}
	    ii:=random(tnodes-1)+1;
	    {writeln('randomnumbertilingvertex ',ii:4,tnodes:5);}
	    {writeln('random ',ii:1);}
	    flst[0]:=0;
	    nuii:=abs(dperp_offset(-nu_offs,nu[ii]));
	    if (((ii<>vpin) and (vpin>=0)) or ((vpin<0) and (ii>tnodes0) or (nuii=1))) then begin
	       ntflp:=ntflp+1;
	       {writeln('aalive ',ii:2);}
	       if(sigmaflip) then begin
		  flipped:=false;
		  if(nv[ii]=2) then check_sigmaflp(ii,lvl,nn_mu,nn_list,flipped);
		  if(flipped) then begin flst[0]:=1; flst[1]:=ii; end;
	       end else begin
		  {writeln('checking-for-flip ',swo:1,ii:3,cpl:2,flipo:2);}
		  check_flp(swo,ii,cpl,flst,flipo,nn_mu,nn_list);
		  if(flst[0]>0) then begin
		     {writeln('dbg-attempt ',swo:1,' ',jj:1,' ',ii:1,' ',flst[0]:1,' ',flst[1]:1);}
		  end;
		  {writeln(swo:1,' ',ii:1,' ',cpl:1,' ',flst[0]:1);}
	       end;
	       {writeln('aaalive');}
	    end;
	 until (flst[0]>0);
      end;
      for o:=1 to flst[0] do fl[flst[o]]:=1; {}
      nuii:=dperp_offset(-nu_offs,nu[ii]);
      de_th:=0.0;
      if (not(deco)) then begin
	 u0:=0.0;
	 for o:=1 to n_obj do obj_st_u[o]:=obj_st[o];
	 echold:=0.0; for o:=1 to n_objch do begin
	    obj_sta_u[o]:=obj_sta[o]; {save}
	    k:=(obj_stn[o]-obj_sta[o]);
	    echold:=echold+obj_ch[o]*k*k;
	 end;
      end;
      if (swo=1) then begin
	 flipped:=(flst[0]>0);
	 bflst[1]:=flipped;
	 if(flipped) then begin
	    nflp:=nflp+1;
	    noflx[flipo]:=noflx[flipo]+1;
	 end;
      end else if (swo=2) then begin
	 flipped:=true; nflp:=nflp+1; bflst[1]:=flipped;
      end else begin
	 {note: sigmaflip occurs in check_sigmaflp before this}
	 if not(sigmaflip) then begin
	    nn_imu:=0;
	    nn_mu:=0; flipped:=true;
	    for j:=1 to flst[0] do begin
	       if flipped then begin
		  if(instances) then begin
		     {we need to reset all instances that contain flipping vertex (vertices)}
		  end;
		  FLIPIT(j,flst,nn_mu,nn_list,flipped);
	       end;
	       bflst[j]:=flipped;
	       if (j=1) then nn_mu0:=nn_mu;
	    end;
	 end;
      end;
      if (tdim<tnodes+nn_mu) then begin
	 writeln('Insufficient array dim for tiling nodes: increase tdim in common.p');
      end;
      if (debug[1] or debug[3] or flip_man) then begin
	 write('attempt ',jj:1,' (node ',ii:1,', nu ',nuii:1,'); flst#=',flst[0]:1,', ');
	 if flip_replay then write('do_flip ',doflip[jj]:1,' ');
	 writeln('flipped: ',flipped);
	 if debug[2] then begin
	    write('nnlist : ',nn_mu:1,' -> ');
	    for i:=1 to nn_mu do write(nn_list[i]:1,' ');
	    writeln;
	 end;
      end;
      if (flipped) then begin
	 {we need 2 node-lists that may differ for multi>1:
	 if io3d=0, the ATOM-lists in the adjacent layers
	 dont need update; but for io3d>1 they do.}
	 nnma:=nn_mu; nn_tmu:=0;
	 if (((lvl>1) or (tlvl>1)) and not(sigmaflip)) then begin {level=0 is no-objective}
	    {note: sigmaflip does this in the main routine}
	    nn0:=nn_mu;
	    if(lvl>0) then add_boundary(lvl-1,nn_mu,nn_list);
	    {if(instances) then begin
	    fff[0]:=0;
	    for i:=nn0+1 to nn_mu do pat_instance_dep(nn_list[i],fff);
	    for i:=1 to fff[0] do begin nn_mu:=nn_mu+1; nn_list[nn_mu]:=fff[i]; end;
	    end;}
	    nodes_union(nn_mu,nn_list);
	    nnma:=nn_mu;
	    if deco then begin
	       for i:=nn0+1 to nnma do COPY_AREC(nn_list[i],tnodes+i);
	    end;
	    if(tlvl>lvl) then begin
	       add_boundary(tlvl-lvl,nn_mu,nn_list);
	       {we need to copy these extra nodes to separate list, since e_upd_buf uses nn_list}
	       for k:=nnma+1 to nn_mu do begin
		  nn_tlist[k-nnma]:=nn_list[k];
	       end;
	       nn_tmu:=nn_mu-nnma;
	       nn_mu:=nnma;
	    end;
	    {for now, nn_list is OK, but will be overwritten after nn_mu in E_update_buff}
	    for i:=nn0+1 to nn_mu+nn_tmu do copy_nrec(nn_list[i],tnodes+i);
	    for i:=nn0+1 to nn_mu+nn_tmu do copy_trec(nn_list[i],tnodes+i);
	    if debug[2] then begin
	       write('ADDED_BOUNDARY ',(nn_mu-nn0):1,'  ');
	       for i:=nn0+1 to nn_mu do write(nn_list[i]:1,' '); writeln;
	    end;
	 end;
	 {for o:=1 to nn_mu do write(nn_list[o]:1,' '); writeln;}
	 {for o:=1 to nn_tmu do write(nn_tlist[o]:1,' '); writeln;}
	 if(instances and (swo<2)) then begin
	    {write('xxx-nnlist ',nn_mu:1,'  ');
	    for k:=1 to nn_mu do write(nn_list[k]:1,":");
	    writeln;}
	    for k:=1 to nn_mu do begin
	       new_pat_instance(nn_list[k],0);
	    end;
	    if(nn_tmu>0) then for k:=1 to nn_tmu do new_pat_instance(nn_tlist[k],0);
	 end;
	 for k:=1 to nn_mu do NEW_NODE_OBJECTS(nn_list[k],1);
	 if(nn_tmu>0) then for k:=1 to nn_tmu do new_node_objects(nn_tlist[k],0);
	 if deco then begin
	    if (swptp=1) then CMPCHL(swo,flst,nn_mu,nn_list);
	    if (swptp=2) then CMPCHLA(swo,flst,nn_mu,nn_list,nla);
	 end;
	 reject:=false;
	 if deco then begin
	    na0:=0; na1:=0;
	    for o:=1 to nnele do dnz[ele[o]]:=0;
	    for i:=1 to n_0 do if (t_0[i]>0) then begin
	       na0:=na0+1;
	       dnz[t_0[i]]:=dnz[t_0[i]]-1;
	    end;
	    for i:=1 to n_1 do if (t_1[i]>0) then begin
	       na1:=na1+1;
	       dnz[t_1[i]]:=dnz[t_1[i]]+1;
	    end;
	    if (swptp=2) then reject:=((nla>0) or reject);
	    if chvar then begin
	       {NOTE: dnz and buflim are array [0..z_ele], while mele_oc and mele_nom are [1..t_ele]}
	       reject:=false;
	       for o:=1 to nnele do
		  if debug[3] then writeln('chvar ',o:1,' ',ele[o]:1,' dnz=',dnz[ele[o]]:1,' ncur=',chcur[ele[o]]:1,' nom=',chnom[ele[o]]:1);
	       for o:=1 to nnele do if not(reject) then begin
		  if(dnz[ele[o]]>0) then begin
		     if(dnz[ele[o]]+chcur[ele[o]]>buflim[ele[o]]+chnom[ele[o]]) then begin reject:=true; {writeln(txtch,'reject_buf+ ',ele[o]:1,' ',dnz[ele[o]]+chcur[ele[o]]:1,' ',buflim[ele[o]]:1)}; end;
		  end else begin
		     if(chcur[ele[o]]+dnz[ele[o]]<-buflim[ele[o]]+chnom[ele[o]]) then begin reject:=true; {writeln(txtch,'reject_buf- ',ele[o]:1,dnz[ele[o]]+chcur[ele[o]]:1,' ',buflim[ele[o]]:1);} end;
		  end;
	       end;
	    end else reject:=(na0<>na1);
	 end;
	 if deco then if ((na0<>na1) or (nla>0)) then begin
	    if(debug[3])then begin
	       writeln('na0 na1 nla ',na0:1,' ',na1:1,' ',nla:1);
	    end;
	    if flip_man then begin
	       if not(chvar) then writeln('WARNING: mismatch in atom numbers. Accepted for debugging.');
	       reject:=false;
	    end else begin
	       if (not(chvar) and reject) then begin
		  nwarn:=nwarn+1;
		  if ((nwarn=1) or debug[1] or debug[2]) then begin
		     writeln('ATOM# MISMATCH (swaptp=',swptp:1,'): flip rejected ',na0:1,' ',na1:1,' ',nla:1,'.');
		     if (nwarn=1) then begin
			export_err:=true;
			writeln('This probably means the decoration file is incomplete or erroneous; else, use "-chp" option to allow chemistry variation');
			writeln('Exporting copas.err, check "fig -fl copas.err"');
			rewrite(txerr,'copas.err');
			writeln(txerr,tnodes:1,' ',multi:1,' 2');
			for i:=1 to 2 do begin
			   for o:=1 to 5 do write(txerr,bai[i,o]:5); writeln(txerr);
			end;
			writeln(txerr,'flipped configuration');
			export_tiling_body_mon(txerr);
			export_params(txerr);
			export_tiling_atoms(txerr);
		     end;
		     if (not(debug[1] or debug[2]) and (nwarn=1)) then writeln('[Will not display this message again]!');
		  end;
	       end;
	    end;
	 end;
	 if reject then begin {undo the flip}
	    {writeln('restore1 ',nn_mu:1);}
	    for k:=1 to tnodes do tvl[k]:=tvls[k];
	    RESTORE_REC(nn_mu,nn_tmu,nn_list,nn_tlist);
	    if export_err then begin
	       {writeln(nwarn:1,' ',export_err);}
	       writeln(txerr,'restored configuration');
	       {writeln('afterneleloop');}
	       if(deco) then begin
		  export_tiling_body_mon(txerr);
		  export_params(txerr);
		  export_tiling_atoms(txerr);
		  close(txerr);
		  export_err:=false;
		  if debug[2] then begin
		     writeln('Current config has been exported into copa.err. [ENTER=continue]');
		     readln;
		  end;
	       end else begin
	       end;
	    end;
	 end;
      end else if not(flipped) then begin {case cpl=1}
	 {flst[0]=2 guarantees we dont do this for sigmaflip}
	 if ((flst[0]=2) and bflst[1] and (swo=0)) then RESTORE_REC(nn_mu0,0,nn_list,nn_list);
      end;
      {writeln(flipped,' ',reject);}
      if (flipped and not(reject)) then begin
	 if (not(deco)) then begin
	    u1:=de_th;
	 end;
	 nuii:=dperp_offset(-nu_offs,nu[ii]);
	 if(swo=0) then begin
	    nflp:=nflp+1;
	    nflx[nuii]:=nflx[nuii]+1;
	 end;
	 if not(deco) then begin
	    echnew:=0.0; for o:=1 to n_objch do begin
	       k:=(obj_stn[o]-obj_sta[o]);
	       echnew:=echnew+obj_ch[o]*k*k;
	    end;
	 end;
	 if (deco) then begin
	    {chpot is common for all species: artificial quadratic term}
	    {DISABLING!}
	    {
	    if ((na0<>na1)) then begin
	     echold:=(natch-natoc)*(natch-natoc)*chpot;
	     natoc:=natoc-na0+na1;
	     echnew:=(natch-natoc)*(natch-natoc)*chpot;
	    end;
	     }
	    dchp:=0;
	    if (chvar) then begin
	       if(flip_man) then write('CHPOT ');
	       for o:=1 to nnele do begin
		  if(flip_man) then write('  ',ele[o]:1,' ',dnz[ele[o]]:1,' ',chp[ele[o]]:1:4);
		  dchp:=dchp+dnz[ele[o]]*chp[ele[o]];
	       end;
	       if flip_man then writeln;
	    end;
	    if fupd then begin
	       if (deco) then begin
		  nn_mu1:=nn_mu;
		  for o:=1 to flst[0] do add2list(flst[o],nn_mu1,nn_list);
		  if (funi=1) then flip_union(feps);
		  E_update_buf(nn_mu,nn_mu1,nn_list,u0,u1);
	       end;
	       if ((debug[3]) or (debug[1])) then begin
		  writeln('upd_fast (u1 u0 dchp du): ',u1:12:6,u0:12:6,dchp:12:6,(u1-u0+dchp):12:6);
		  dux:=u1-u0+dchp;
		  if (debug[3]) then begin
		     E_update_thru_nodes(nn_mu,nn_mu1,nn_list,u0,u1);
		     writeln('upd_all  (u1 u0 dchp du): ',u1:12:6,u0:12:6,dchp:12:6,(u1-u0+dchp):12:6);
		  end;
	       end;
	       u1x:=u1; u0x:=u0;
	    end;
	    {dramax:=0.1;
	    if (relmcs>0) then begin
	    for i:=1 to n_pxyz do
	    for j:=1 to nn_pxyz[i] do dpxyz0[i,j]:=dpxyz[i,j];
	    relax_para(relmcs,xgamma,dramax,u1);
	    end;}
	 end;
	 if flip_replay then begin
	    if not(skip_reject and (doflip[jj]=0)) then export_replay(jj,ii);
	 end;
	 if flip_rec then write(txfout,ii:1,' ',(tpot+u1-u0+dchp):1:6,' ',(u1-u0+dchp):1:6);
	 du2:=(u1-u0+dchp)*(u1-u0+dchp);
	 du0:=beta[kk]*(u1-u0+echnew-echold+dchp);
	 if flip_man then begin
	    dumt:=dumt+u1-u0+dchp;
	    writeln('FLIP_MAN: dU=',(u1-u0+dchp):1:5,' [eV]');
	 end;
	 force_reject:=(flip_replay and (doflip[jj]=0));
	 force_accept:=((flip_replay and (doflip[jj]=1)) or flip_man);
	 if (flip_replay and (jj=til0)) then begin
	    writeln('tiling ',jj:1,' -> tiling.rec');
	    for i:=1 to tnodes do fl[i]:=0;
	    fl[ii]:=jj;
	    export_tiling_mon('tiling.rec');
	 end;
	 if (flip_replay and (jj>=til1) and (jj<=til2)) then begin
	    for i:=1 to tnodes do fl[i]:=0; fl[ii]:=jj;
	    writeln('tiling ',jj:1,' -> tilings.rec');
	    export_tiling_body_mon(txti);
	 end;
	 try_reject:=(not(force_accept) and ((du0>eeps) or force_reject));
	 if try_reject then begin
	    if(du0>20) then x:=0 else x:=exp(-du0);
	    y:=randu(seed);
	    if ((y>x) or force_reject) then begin {reject flip}
	       nflp:=nflp-1;
	       if(swo=0) then begin
		  nflx[nuii]:=nflx[nuii]-1;
	       end;
	       if(swo=1) then begin  {old object flip}
		  noflx[flipo]:=noflx[flipo]-1;
	       end;
	       for o:=1 to n_objch do obj_sta[o]:=obj_sta_u[o];
	       for o:=1 to n_obj do obj_st[o]:=obj_st_u[o];
	       {writeln('restore2 ',nn_mu:1);
	       for o:=1 to nn_mu do write(nn_list[o]:1,' '); writeln;}
	       for k:=1 to tnodes do tvl[k]:=tvls[k];
	       RESTORE_REC(nn_mu,nn_tmu,nn_list,nn_tlist);
	       if (debug[1]) then write('dU=',(u1-u0):1:6,' dchp=',dchp:1:6,' dUtot=',(u1-u0+dchp):1:6,' : REJECT');
	       if flip_rec then writeln(txfout,'  0');
	       if (deco) then begin
		  natoc:=natoc+na0-na1;
		  {for i:=1 to n_pxyz do
		  for j:=1 to nn_pxyz[i] do dpxyz[i,j]:=dpxyz0[i,j];
		  if (relmcs>0) then begin
		  ipos:=1; iposd:=1;
		  for i:=1 to tnodes do
		  for j:=1 to nan[i] do
		  para2xyz(n_para[i,j],n_a[i,j,1],ipos,iposd,n_r[i,j]);
		  end;}
	       end;
	    end else begin {accept}
	       for o:=1 to nnele do chcur[ele[o]]:=chcur[ele[o]]+dnz[ele[o]];
	       if flip_rec then writeln(txfout,'  1');
	       if (debug[3] and (deco)) then acer:=acer+u1-u0+dchp-dux;
	       tpot:=tpot+u1-u0+dchp;
	       mystat(kk,lab,tpot);
	       if ((jj mod save_e)=0) then begin
		  write(recu,(u1-u0):1:6,' ',kk:1,' 4 nat=',natoc:1,' nui=',nuii:1,' ');
		  if(chvar) then begin
		     for o:=1 to nnele do begin
			if(o>1) then write(recu,',');
			write(recu,ele[o]:1,":",dnz[ele[o]]:1);
		     end;
		     {write(recu,' dchp=',dchp:1:5,' ');}
		  end;
		  writeln(recu);
	       end;
	       if (debug[1]) then begin
		  write('dU=',(u1-u0):1:6);
		  if(chvar) then write(' dchp=',dchp:1:5);
		  write(' : ACCEPT+');
	       end;
	    end;
	 end else begin {accept}
	    for o:=1 to nnele do chcur[ele[o]]:=chcur[ele[o]]+dnz[ele[o]];
	    if flip_rec then writeln(txfout,'  1');
	    if (debug[3] and (deco)) then acer:=acer+u1-u0+dchp-dux;
	    tpot:=tpot+u1-u0+dchp;
	    MYSTAT(kk,lab,tpot);
	    if ((jj mod save_e)=0) then begin
	       write(recu,(u1-u0):1:6,' ',kk:1,' 4 ',natoc:1,' ',nuii:1);
	       if(chvar) then writeln(recu,' dchp=',dchp:1:5) else writeln(recu);
	    end;
	    if (debug[1]) then
	       write('dU=',(u1-u0+dchp):1:6,': ACCEPT-');
	 end;
	 if debug[1] then writeln(' : U=',tpot:1:5);
      end;
      {if (deco) then x:=tpot/natoc else x:=tpot/tnodes;}
      if (deco) then x:=tpot else x:=tpot/tnodes;
      if (x<emin) then export_config_if(emin,x,beta[kk],sname[lab]);
      if (debug[3] and (deco)) then
	 writeln('Accumulated error: ',acer:1:7);
      {writeln('FFF ',emin:1:4,' ',x:1:4);}
      for o:=1 to flst[0] do fl[flst[o]]:=0;
      {writeln('gggaaa');}
   end; {loop flips}
   rate:=nflp/ntflp;
   if (flip_man) and (nseq>1) then begin
      writeln('FLIP_MAN: dU(TOT)=',dumt:1:5);
   end;
  { if(flip_man) then begin
      writeln('copa.out is restored configuration, copa.tmp flipped.');
      rewrite(txerr,'copa.tmp');
      writeln(txerr,tnodes:1,' ',multi:1,' 2');
      for i:=1 to 2 do begin
	 for o:=1 to 5 do write(txerr,bai[i,o]:5); writeln(txerr);
      end;
      writeln(txerr,'restored configuration');
      export_tiling_body_mon(txerr);
      export_params(txerr);
      export_tiling_atoms(txerr);
      close(txerr);
      RESTORE_REC(nn_mu,nn_tmu,nn_list,nn_tlist);
   end;}
end; { mc_tiling }

end.





